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0.  Introduction 


The  approximation  of  partial  differential  equations  by  the  finite  element  or 
finite  difference  methods  often  leads  to  large  sparse  linear  systems. 
Traditionally,  one  employs  basic  iterative  methods  to  solve  them.  However, 
basic  iterative  methods  exhibits  slow  convergence  if  the  meshsize  is  small. 
This  is  due  to  the  fact  that  basic  iterative  methods  can  not  damp  out  low 
frequency  modes  of  the  errors.  To  remedy  this  disadvantage,  multigrid  methods 
combine  basic  iterative  methods  with  other  methods  that  are  complementary. 
One  of  the  reasons  that  accounts  for  the  effectiveness  of  multigrid  methods 
seems  to  be  the  idea  of  approximating  the  solution  of  a  large  system  from  a 
subspace  whose  dimension  is  small.  Finding  such  a  subspace  is  by  no  means  an 
easy  task.  In  this  report  we  investigate  the  residual  based  reduced  basis  method 
(RRBM)  for  solving  large  sparse  symmetric  positive  definite  systems.  In  a 
sense,  it  belongs  to  the  class  of  two-grid  methods,  although  no  geometric  grids 
are  involved.  The  outline  of  this  report  is  as  follows.  Section  1  is  devoted  to 
the  convergence  behavior  of  general  projection  processes.  A  generic  bound  is 
provided  for  a  class  of  pseudoresidual-based  projection  processes.  In  Section  2, 
we  give  mild  conditions  that  ensure  the  convergence  of  general  additive 
correction  procedures.  In  Section  3,  we  describe  the  RRBM  and  prove  a 
convergence  theorem  for  it.  Section  4  is  devoted  to  the  practical 
implementation  of  the  RRBM.  In  Section  5  we  relate  the  RRBM  to  the 
preconditioned  conjugate  gradient  method  under  a  certain  condition.  It  is  shown 
that  the  RRBM  is  more  flexible  than  the  preconditioned  conjugate  gradient 
method.  In  Section  6  we  extend  some  of  the  previous  results  to  the  case  of 
nonsymmetric  linear  systems  having  positive  definite  symmetric  part.  A  rather 
general  convergence  analysis  is  given  there.  Finally,  Section  7  is  devoted  to  the 
weight  selection  principle  for  the  linear  convective  equation  in  Rn.  The 
principle  is  based  on  the  ability  of  hybrid  difference  methods  to  conserve  th^ 
discrete  weighted  energy. 
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1 . 

Additive  Correction  Scheme 

j 

Let  A 

be  a  symmetric  positive  definite  (SPD) 

matrix 

of  order  n. 

Consider  the 

linear  system  of  equations 

Au  =  b. 

(1.1) 

The 

general 

additive  correction  scheme  for  solving  (1.1) 

can  be 

described  as 

follows. 

ALGORITHM  1:  Let  u(°)  be  aiven.  Choose  a  fixed  inteaer  u.  For 

k  = 

1,2 .  do 

step  1. 

Initialization; 

w(0)  s  u0c) 
k 

Step  2. 

Presmoothing: 

w(j)  =  Gw,(j'1}  ♦  0" 1  b,  j  =  1.2 . u. 

k  k 

(1.2a) 

Step  3. 

Defect  Computation: 

d*°*  =  b  -  Awfu^  =  A(u  -  w<u))  s  Ae[u'. 
k  k  k  k 

(1.2b) 

Step  4. 

Additive  Correction: 

i/k+1  ^  =  w^u^  ♦  E 
k  R 

(1.2c) 

where  e  is  a  computationally  inexpensive  approximation  of  e^  . 

K  K 


Some  remarks  are  in  order.  Steps  1  and  2  in  the  above  algorithm  consists 
of  the  usual  basic  iterative  method  with  the  iteration  matrix  G  and  the 
splitting  matrix  0.  It  is  well  known  that  basic  iterative  methods  applied  to 
the  linear  systems  arising  from  discretization  of  partial  differential  equations 
are  not  effective  once  the  high  frequency  modes  have  been  damped  out.  Hence  one 
stops  after  a  certain  number  of  steps  of  the  basic  iterative  method.  At 

this  juncture,  tne  error  will  lie  mainly  in  the  subspace  of  low 

frequency  modes.  Intuitively,  this  means  the  error  can  be  well 

approximated  oy  solving  an  approximating  system  of  the  linear  system  (1.2b)  in 
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a  smaller  subspace.  The  solution  e  (R  for  'reduced')  is  then  added  to 

The  success  of  Algorithm  1  depends  whether  one  can  judiciously  choose  the 
above-mentioned  subspace.  In  this  report  we  confine  ourselves  to  projection 

processes  in  which  e  is  the  projection  of  onto  a  certain  subspace. 

H  K 

We  now  study  the  effects  of  steps  3  and  4. 

Let  w  be  an  approximation  of  the  solution  u  to  the  system  (1.1).  Let 

S  =  [s i  ,s2 sm]  be  a  full  rank  matrix  consisting  of  column  vectors  sj, 

1  <_  i  <  m.  The  following  describes  how  to  get  a  new  approximation  u  from 
w  . 


Galerk  in 


lation 


Id  s  b  -  Aw  =  Au  -  Aw  a  Ae.  (1.3a) 

Find  e  e  Rg(S)  =  the  column  space  of  S  such  that 

R 

ST(d  -  Ae  )  =  0.  (1.3b) 

H 

U  =  W  ♦  E  .  (1.3c) 

H 

The  above  is  equivalent  to 

SR  =  S(STAS)'1STAe.  (1.4a) 


e  =  u-  u=  e-e  =e  -  S(STAS)_1  STAe.  (1.4b) 

H 

Let  (x,y)g  =  xTAy.  the  energy  product  associated  with  the  matrix  A.  It  is 

easy  to  see  that  P$  s  S(STAS)*1STA  is  the  orthogonal  projector  onto  Rg(S). 

Let  ii •  ii ^  be  the  energy  norm  induced  by  (•,•)£•  Then  from  (1.4b)  we  have 

lien  =  ii(I  -  P  )e H  =  min  lie  -  zii  (1.5a) 

E  s  E  ZeRgfS)  E 

or 


iiu  -  uii  =  min  iiv  -  uii.  (1 .5b) 

E  vewRg(S) 

It  is  known  that  if  Rg(S)  =  m ) (d.A)  s  Span{d.Ad A^'^dl,  then  (1.5a) 

results  in  a  bound  involving  a  Chebyshev  polynomial  of  the  first  kind.  That  is 
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where  Tm  is  the  Chebyshev  polynomial  of  degree  m  and  [a,bj  is  the  smallest 
interval  containing  the  spectrum  of  the  matrix  A.  Note  that  the  conjugate 
gradient  method  falls  into  this  category.  We  prove  below  t ha L  a  generic 
inequality  exists  for  all  projection  processes  (1.4)  if  d  e  Rg(S). 


LEMMA  1.1  [Lemma  3.1,  D.  Braess,  1981]  Let  the  Hilbert  space  U  be  a 
direct  sum  of  its  subspaces  V  and  W.  Assume  that  there  is  a  Z  <  1  such 
that  a  strengthened  Cauchy  inequality  holds: 

|  (v,w)  |  <  Z\\  VIIIIWII.  V  e  V.  W  e  W. 

If  Py/u  =  0,  then 

IIU  -  Pvull  <  yilUII. 

Here  Pw  and  Py  are  orthogonal  projectors  to  W  and  V,  respectively. 


THEOREM  1.1.  Let  e  and  e  be  as  in  (1.4).  Assume  that  the  defect 
d  e  Rg(S).  Then 

■ .  /  *  \  « 

(1.6) 


lie,!  <  MA)  -  1. 
E  “  k(A)  ♦  1 


E’ 


where  k(A)  =  it  Ah  2  II  A- 1 II 2  and  imi2  b  the  2-norm  of  A. 

Proof.  We  set  various  spaces  corresponding  to  Lemma  1.1  as  follows. 
U  =  Rn  endowed  with  the  energy  product. 

V  =  Rg(S),  W  =  Rg1(S)  =  the  orthogonal  complement 
of  V  with  respect  to  the  Euclidean  inner  product. 


Obviously, 


Rn  =  V  ©  W. 


Furthermore,  Py/e  =  0  since  Ae  =  d.  Also,  Py  =  S(STAS)_  1 STA.  The  constant 
Z  in  Lemma  1.1  can  be  determined  as  follows.  For  v  e  V  and  w  e  W. 
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(v,w)£  -  vtAw 


=  VTAW  -  o<VT W 


vt(A  -  od)w 


V<*  >  0 


=  (A1  /2v)T  CA-1/2(A  -  o<i)A-1/2)A1/2w. 


Hence 


]  (v.w)e  |  <  II  VII E  II w II £  III  -  cXA'1!^. 

Let  0  <  X,  £  X2  <  ...  <  Xn  be  the  eigenvalues  of  A.  Now 

7S  =  min  ii  I  -  cxA'1  ii 

o<>0  2 

=  min  max  |l  -  <*X  | .  o(A_1)  =  the  spectrum  of  A' 
°<>0  Xe<J(A*’) 

=  min  { j  1  -  cxX*1  | . |  1  -  cxX"1  |  } 

o<>0  1  n 

X  -  X, 
n  1 

X  ♦  X 
n  1 

=  Ck(A)  -  l)/(k(A)  ♦  i). 

The  conclusion  of  the  theorem  then  follows  easily  by  noting  that 

e  =  (I  -  Pv)e. 


We  remark  that  Theorem  1.1  covers  the  steepest  descent  method  and  the 
conjugate  gradient  method.  Inequality  (1.6)  for  the  steepest  descent  method  is 
usually  derived  from  Kantorovich  inequality  (see,  e.g.  Luenberger,  1984). 

Theorem  1.1  can  be  generalized  as  follows. 


THEOREM  1.2.  Let  e  and  e  be  as  in  (1.4).  Assume  that  d  e  Rg(S).  Let 
C  be  a  SPD  matrix  such  that  C(Rg(S))  c  Rg(S).  Then 

__  L 1  /2  a  /■**  1  /2\  i 


m eii  <  - QV  — u  lien,.. 

E  k(C'1/2AC'1/2)  ♦  1  E 

Proof.  The  corresponding  subspaces  are  the  same  as  in  Theorem  1.1.  The 
only  detail  changed  is  the  choice  of  7S.  Let  v  e  V  and  w  e  W.  Then 
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(v,W)E  =  wtAv 
=  wtAv 


<*WTCV 

=  WT(A  -  <*C). 


Vo<  >  0 


As  before,  we  have 

|(v.w)e  |  <  »i w H £  ii vii £  ill  -  ofA-1  ^CA'1 /2ii2. 

Set  B'1  =  A-1 /2CA"1 /2.  For  a  SPD  matrix  D  we  denote  the  eigenvalues  of  D 
as  0  <  X^D)  <  X2(D)  <  ...  ±  Xn(D).  Define  the  constant  Z  in  Lemma  1.1  as 

..  >tm-i  VB)-Xim 


But 


and 


k(B)  *  1 


X  (B)  ♦  X  (B) 
n  l 


X  (B)  =  p(B)  =  the  spectral  radius  of 
--  p(A1/2C-]A,/2) 

=  p(C_1A) 

=  p(C~1/2C"1/2A) 

=  p(C'1/2AC'1/2), 


Xt(B)  =  1/(1/X  (B)) 

=  1 /p(B” 1 ) 

=  1/p(A',/2CA~,/2) 

=  l/pfCA'1) 

-  1/p(C1/2A'1C1/2) 

=  1/p((C'1/2AC'1/2)-1) 

--  x,(c1/2ac-1/2). 


QED. 


IHiOREM-  J  t3,. 


Let 


S  = 
lien 


Cd  e  Rn* '  and  0  x  d. 

<  KiQl^hQ.] lle 

E  k(C1/2AC1/2)  ♦  1 


If  c 

E' 


Proof.  In  reference  to  Lemma  1.1,  we  set 

V  =  SpanfCdh 


is  SPD  then 
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W  =  {w  |  wTd  =  0} 
U  =  Rn. 


It  ts  easy  to  see  that 
and 


U  =  V  ©  W, 


Py^e  =  0. 

Let  v  g  V  and  w  e  W.  v  e  V  implies  that  v  =  0Cd  for  some  0.  Thus 


(v,w)  =  vtAw 

=  vTAw  -  c*0dTw.  V<*  >  0 

=  vTAw  -  o^ttT1  )Tw 
=  vT(A  -  olC'1  )w. 

Now  the  rest  of  proof  is  just  as  in  Theorem  1.2,  we  omit  the  details. 

OED. 


If  we  take  C  =  1,  the  identity  matrix,  we  recover  the  steepest  descent 
method.  If  a  has  Property  A.  then  we  can  take  C  =  D"1,  where  D  is  the 
diagonal  part  of  the  matrix  A.  It  is  well  known  in  this  case  that 
D-1/2AD-1/2  has  smaller  condition  (see,  e.g.  Young  1971,  p.  214).  In  general, 
we  can  take  C  =  Q‘1,  where  0  is  the  splitting  matrix  of  a  basic  iterative 
method  applied  to  the  system  (1.1). 


2.  Convergence  Analysis 


We  now  turn  to  convergence  analysis  of  Algorithm  1.  Since  we  are 
interested  in  the  relationship  between  two  consecutive  iterates,  we  shall  drop 
the  subindex  k  appearing  in  (1.2a-1.2c).  One  cycle  of  Algorithm  1  with 
Galerkin  approximation  is  as  follows. 

Given  u^. 

1.  Initialization; 

yy(0)  =  u(k). 


W  ^  =  GW^~^  ^  ♦ 


2.  Presmoothing: 


O' 1  b,  j  =  1.2 . u. 


(2.1a) 


9 


3.  Defect  Computation: 


d  =  b  -  Aw^  =  A(u  -  w^uh  =  Ae^u\ 

(2.1 b) 

4.  Additive  Correction: 

Let  a  full 

rank  matrix  S  e  Rn*m  be  given. 

(4. a)  Solve 

StASz  =  STd. 

(2.2) 

(4.b)  Set 

£  =  SZ. 

R 

(2.3) 

(4.c)  Set 

=  W(u)  *  £  . 

R 

(2.4) 

Define  the  k-th  error 

vector  as 

e^  =  u  -  u<k). 

(2.5) 

It  is  easy  to  see  that 

e(u)  =  Gue(k) 

(2.6) 

and 

e(k*1}  =  £(u)  -  £  . 

R 

(2.7) 

From  (1 ,4b)  with  e  = 

e^*1)  and  e  =  we  have 

e(k  +  1)  =  (l  -  S(StAS)'1StA)e(u) 

(2.8) 

=  (I  -  S(STAS)_1STA)Gue(k). 

(2.9) 

By  (2.9), 

A 1  /2e(k  ♦  1 ) 

=  (I  -  A1/2S(STAS)-1sTA1/2)A1/2GuA-1/2A1/2e(k). 

Note  that 

P,  =  A1/2S(STAS)'1STA1/2 

(2.10) 

is  the  orthogonal  projector  (with  respect  to  the  Euclidean  inner  product)  onto 
the  range  of  A1/2S.  Hence 

A1/2e(k*1}  =  (i  -  Pl)A1/2GuA*1/2A1/2e(k).  (2.11) 


whence 


1  0 


lle(k  +  1)llE  <  ill  -P^i  /iGuli£  »e(k)liE 

=  ll Gu II  lle{k)llE.  (2. !  2) 

Recalling  that  G  =  I  -  Q'^A  for  a  basic  iteration  matrix  G,  we  immediately 
have  the  following  conclusion. 

THEOREM  2.1.  If  0  is  symmetric  and  p(G)  <  1.  then  any  sequence 
produced  by  Algorithm  1  with  Galerkin  approximation  converges  to  the  true 
solution  u. 

Proof.  By  (2.12), 

iie(k* 1  Jn E  <  ii Guii ^  ne(k)nE 

<  llGllj?  lie(k)llr  (2.13) 

£  £ 

-  ni  -  Q~ 1  An E  iie^ii 
=  ill  -  A,/20',A1/2ii  “  ne(k)np 
=  II  A1  /2GA'  1  /2|!  2  lie(k)llE 
=  pu(G)  ne^k)liE. 

QED. 

THEOREM  2.2.  If  we  use  the  Gauss-Seidel  method  for  the  presmoother  in 
Algorithm  1  with  Galerkin  approximation,  then  any  secuence  produced  by 
Algorithm  1  is  convergent  to  the  true  solution  u. 

Proof.  According  to  a  theorem  in  Young  [1971,  p.  79], 

ii  G  ll  £  <  1  if  and  only  if  0T  ♦  0  -  A  is  SPD. 

Noting  (2.13)  and  applying  the  theorem  to  the  Gauss-Seidel  iteration  matrix,  we 
obtain  Theorem  2.2. 

OED. 


Now  by  (1.5a)  with  e  =  and  e  =  we  have 


Vz  e  Rg(S). 


1  1 


min  ii 

ZeRg(S) 


e(u)  -  2ii. 


lGue(l°  -  zil. 


(2.14) 


Inequality  (2.M)  suggests  that  we  choose  Rg(S)  to  be  the  subspace  spanned  by 
{Gle^h  This  is  what  we  shall  do  in  the  next  section. 


3.  Reduced  Basis  Technique 

From  (2.9)  and  Section  1.  we  know  that  e^*1)  is  the  orthogonal 
projection  of  Gue^^  onto  the  Rg(S).  Hence  the  column  vectors  of  the  matrix 
S  can  be  replaced  by  any  basis  of  the  range  Rg(S).  What  matters  is  the 
subspace  spanned  by  the  column  vectors  of  S.  It  is  now  clear  that  once  we 
specify  a  subspace  of  Rn.  Algorithm  1  is  then  completely  defined. 

From  now  on  we  shall  denote  one  cycle  of  Algorithm  1  with  Galerkin 
approximation  by  RBM(G.Sr.u).  Here  RBM  stands  for  the  reduced  basis  method 
and  Sr  is  the  reduced  subspace.  Consider  the  following  choice  of  Sr,  which 
was  first  proposed  by  Porsching  [1990], 

Pseudoresidual  Based  RBMCG.Sp ,u): 


Sr  =  span{8^  ).8^2) . 8^u)}_ 

where 

J(i)  :  w0)  .  u(k)i  j  =  |  2 . u.  (3.1  ) 

Let  k(u)(x.B)  a  span{x.Bx . Bu_1xh  x  e  Rn.  Be  R1"1*1"1. 

LEMMA  3.1.  Define  5^)  =  Gu^)  ♦  0_1b  -  u^k).  Suppose  that  i  -  G  is 
invertible,  then 

S  =  Spanf^15 . S<u)} 

H 

=  K(u)(S(k),G)  =  K(u)(S(k)J  -  G). 


Proof.  rW  =  w(i}  -  u  *  u  -  u 


=  -G^e^  ♦  e^ 

=  (I  -  Gj)e(k). 

s(k)  _  Qu(k)  t  Q-1b  .  u(k) 


=  Gu^k  ^  ♦  0_1Au  -  u^ 


(3.2) 


Hence 


(I  -  G)e( 


S<j)  =  (1  -  Gl)(l  -  G)' 1  R^ 


:  (I  -  G ) ( I  ♦  ...  ♦  GJ " 1  )(I  -  G)' 1  S 
=  U  ♦  G  ♦  ♦  Gi-1  )Stki. 


-l £(k) 


(3.3) 


we  remark  that  the  condition  1  -  G  is  invertible  holds  for  any  completely 
consistent  basic  iterative  method.  This  is  true,  for  most  commonly  known  basic 
iterative  methods.  6^)  is  often  called  pseudoresidual  vector.  From  now  on, 
we  shall  denote  the  pseudoresidual-based  RBM(G.Sr.u)  as  RRBM(G.Sr.u). 


THEOREM  3.1  [Chou  and  Porsching,  1988]  Given  a  RRBM(G.Sr.u)  with 
1  -  O' 'a  defining  a  completely  consistent  presmoother,  if  Q  is  SPD  then 

a)  G  has  real  eigenvalues  { ja j }.  j  =  1,2 . n,  which  may  be  ordered  as 

Pi  <  P2  <  -  <  PN  <  1- 

b)  ii k * 1  )|i£  <  T\(u)ne^)n £,  where 

■q(u)  =  22fu/2/(1  ♦  2fu)  <  1. 


(1  -  V  1  -  a2  )/(l  ♦  v  1  -  o2  ). 


0  =  (^-^)/(2-^-m,). 


(3.4) 


Proof.  Statement  (a)  follows  easily  by  noting  that 

Al/2GA-»/2  =  ,-A1/2Q-1A1/2. 


Statement  (b)  can  be  proved  as  follows.  By  (2.14)  with  Rg(S)  =  Sr, 

u-1 

ne(k  +  1)llc  <  iiGue(k)  -  X  eX.Gj5(k)l! 

E~  U  1  E 

u-1 

=  llGue(k)  -  ^  erf.Gj(l  -  G)e(k)ll E 
j=0  ] 

=  np(G)e^n  .  for  all  polynomials  p(x)  of  deg  <  v  and  p(l)  =  1. 
<  llA1/2p(G)A_1  /2ll2lie(k)llE 
=  p(p(A1/2GA‘1/2))lle(k)llE 
=  p(p(G))iie(k)ll  E- 

Hence  ii e^k ■" 1  ^ n  <  min  max  |p(X)|ne^ii  .  Define 
E  deg  p  <  u  p,<X<p  E 

P(l )  =  l 

r[(u)  =  min  max  |p(X)|.  The  conclusion  of  the  theorem  now  follows 
deg  p  <u  p  <X<u 
p(l)=l  N 

from  the  well-known  Chebyshev  minimax  theorem  [Young,  1971,  p.  302], 

QED. 


We  now  give  a  theorem  that  characterizes  u(k+1)  knowing  u(k). 

THEOREM  3.2.  Given  a  RRBM(G.Sr.u)  with  initial  iterate  u(k).  Then 
u{k>1)  minimizes  the  quadratic  functional  F(w)  a  iiw  -  uii2  over  the  linear 
manifold  u^k^  ♦  Sr.  That  is 

llu(k^  -  u»E  <  IIW  -  UII E  Vw  e  u^  *•  K^u^(S(k\l  -  G).  (3.5) 

Proof.  By  (2.14). 
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ii e^k * 1  ^ii  =  min  nGue^  -  211 

£  “SR  E 

=  min  iie(u)  -  211 

ZeS„  E 


=  mm  11  u  -  w 

Ze  S„ 


(u) 


211. 


mm  11  w  -  uii. 


WeW^Sp 

The  theorem  will  be  proved  if  we  can  show  that 


w(u)  +  S  =  uEk)  ♦  S  . 
R  R 


Now 

v 

w(u)  =  ^  (w^}  -  w(j'n)  *  w(0) 

=  ^  (Gw(i_1)  ♦  Q"1  b  -  w^'15)  ♦  u(k) 

ji,1 

=  ^  (1  -  G)(u  -  w(j-n)  ♦  u<k> 

=  X  (I  '  G)Gj'1e(k)  u(k) 

V 

=  ^  (I  -  G)G^“ 1  (I  -  G)_1S(k)  ♦  u(k) 
j  =  1 


=  ^jT  G j ' 1  5<k)  ♦  u(k). 
j  =  1 

Hence  w(°)  -  u(k)  e  Sr. 


(3.6) 


(by  2.1a) 


(3.7) 

QED. 


COROLLARY  3.2.1 .  RRBM(G.Sr.u)  with  an  initial  iterate  u(k)  is 
equivalent  to  v  steps  of  conjugate  gradient  method  with  initial  iterate  u(k). 
provided  that  G  =  I  -  A. 

Proof.  Note  that 


By  (3.5). 


S  =  Span{S(k3.(!  -  G)8(k) . (I  -  G)u' 1  8(k)f 

H 

=  Span{8(k),AS(k) Au'18(k)f 

=  Span{r^k3.Ar^k3 . Au'Vk3},  r^k3  =  b  -  Au^k3. 

nu(k  +  l3  -  uiig  <  ii w  -  ui> £  Vw  e  u(k3  ♦  K^u)(r^k3,A). 

OED. 


Corollary  3.2.1  suggests  a  procedure  for  implementing  RRBM(G.Sr.u)  based 
on  Theorem  3.2.  Recall  that  minimization  of  a  quadratic  functional  induced  by  a 
SPD  matrix  B  over  a  k-plane  can  be  achieved  by  minimizations  of  the  quadratic 
functional  along  k  B-conjugate  directions  on  the  k-plane  (see  Hestenes,  1980,  p. 
101).  Thus  we  can  implement  RRBM(G.Sr.u)  by  the  conjugate  direction  method 
and  avoid  explicit  implementation  of  (2.1a).  (2.1b).  (2.2)-(2.4)  altogether. 


4.  Implementation  of  RRBM(G.Sr.u) 


Perform  the  following  steps  for  n  =  1 . v  -  1. 

W  (0  )  s  y(k  ) . 

w(n+ 1 )  =  w(n)  *  xnp(h). 

_(n)  .  /5"”  if  "  =  0 

|s(n)  ♦  o^"'13.  n  =  1.2 . V  -  1, 

<*n  =  -(8(n).AP(n-1))/(p(n-1)>Ap(n-1)). 

s(n)  =  Gw^n3  ♦  Q-1  b  -  w^n3. 

Xn  =  (pCnJ.osCnJj/fpfnJApfn)), 

=  (ptn),r(n3)/(p(n3,Ap(n3), 
r(n)  =  b  -  Aw(n3, 


G  =  I  -  Q"1 A.  (4.1) 

Set  u^k  * 33  =  w<u3. 

Here  (-.-)  denote  the  Euclidean  inner  product.  <*n  is  obtained  by  insisting  that 
(p('3,Ap03)  =  o.  i  *  j.  The  Xn  is  obtained  as  follows. 


and 


Since  w<n*1)  =  w'r')  *  Xnp(n5,  it  is  not  hard  to  see  that  , 

S<n*1)  ♦  \np(n)  =  s(n)  ♦  xnGp(n) 

r(nt,)  -  r!n)  --  -X  Ap(nl.  (4.2) 

n 

Furthermore,  w(n+1)  is  the  minimizer  along  the  direction  p(n).  Hence 
(r(n*  1  )  p(n))  =  o.  This  together  with  (4.2)  implies 

Xn  =  (r(n).p(n))/(p(n),Ap(n)). 

5.  Relationship  between  RRBM(G.SRlu).  G  =  1  -  0-1A  and 
Preconditioned  Conjugate  Procedures 

Let  (I  -  G)u  =  c.  where  c  =  0" 1  b.  Furthermore  let  I  -  G  be 
symmetrizable  with  a  symmetrization  matrix  W  [Hegaman  and  Young,  1981]. 
Define 

4  ■  W(l  -  G)W'l , 
u  =  Wu, 
b  =  Wc. 

Then  Au  =  b  is  called  the  preconditioned  system  with  respect  to  W.  One  can 
apply  conjugate  gradient  method  to  this  system.  Then  the  following  formulae 
are  obtained  [Hegeman  and  Young,  1981,  p.  146]. 

u<0>  is  arbitrary 
u(n+D  =  u(n)  ♦  Xnp^n),  n  =  0.1 . 

(n)  /s(0) •  n  =  °- 

P  '  \6(n)  ♦  c*  p(n_1),  n  =  1.2 . 

nr 

o<n  =  -(WsM.WO  -  G)p(n-n)/(wp(n-1).W(!  -  G)p<n-D).  n  =  1.2 . 

Xn  =  (Wp(n).W5(n))/(WpW,W(I  -  G)p(n>),  n  =  0,1,2 . 

Here 

s(n)  _  Gu(n)  +  c  .  u(n)  (5 . 1 ) 


THEOREM  5.1.  Let  RRBM(G,SR,u)  be  such  that  G  =  I  -  O’1  A  has  a  SPD 
splitting  matrix  0.  Then  RRBM(G,SR,u)  with  an  initial  iterate  u  is 


equivalent  to  the  preconditioned  conjugate  gradient  procedure  (5.1)  .with 
W  =  Ql  /2,  u(0)  =  u  and  the  upper  limit  of  n  =  v.  In  particular, 

RRBM(G,Sr.oo)  is  equivalent  to  the  preconditioned  conjugate  gradient  procedure, 
provided  that  the  same  initial  iterate  is  taken. 

Proof.  Compare  (5.1)  with  (4.1).  (See  Luenberger  [1984.  p.  245].) 

QED. 

We  remark  that  if  Q  is  not  SPD.  the  two  methods  are  not  equivalent.  On 
the  other  hand.  (4.1)  is  still  applicable  for  nonsymmetrizable  cases.  For 

instance.  Theorem  2.2  guarantees  the  convergence  of  U  RRBMtG.S  u)  when 

v  R 

G  is  the  iteration  matrix  of  the  Gauss-Seidel  smoother. 

6.  Contraction  Numbers  in  the  Positive  Real  Case 

In  this  section  we  derive  generic  contraction  numbers  for  a  class  of 
additive  correction  methods  based  on  orthogonal  projection.  The  only  assumption 
on  the  range  of  the  projector  is  that  it  contains  the  residual.  This  generalizes 
the  results  of  Section  1  (SPD  case).  However,  the  contraction  number  obtained 
here  is  not  as  sharp  as  the  earlier  one.  The  general  approach  below  provides 
convergence  results  for  restarted  generalized  conjugate  gradient  methods  under  a 
variety  of  conditions. 

Consider  the  problem  of  determining  u  e  Rn  such  that 

Au  =  f,  (6.1) 

where  A  e  Rnxn  is  invertible,  and  f  e  Rn. 

Certain  popular  iterative  methods  for  the  solution  of  (6.1).  such  as 
generalized  conjugate  gradient  (GCG)  methods  (Hageman  and  Young  [1981.  p.  339], 
Elman  [1982],  Saad  and  Schultz  [1985],  Vatsya  [1988])  and  multigrid  methods 
(Hackbusch  [1985],  McCormick  [1987])  incorporate  into  their  overall  solution 
strategies  an  additive  correction  algorithm  of  the  following  type:  Given  an 
approximation  u0  of  u  and  an  m  dimensional  subspace  Sm  c  Rn: 

1 .  Compute  the  residual 


r  =  f  -  Au  . 
0  0 


(6.2) 


) 
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2.  Observe  that  the  error  e0  =  u  -  u0  satisfies 

Aeo  =  r0.  (6.3) 

3.  Compute  an  approximation  eQ  e  of  eQ. 

4.  Form  a  corrected  approximation  of  u. 

u,  *  u0  *  V  (6'4) 

The  key  to  this  additive  correction  phase  is.  as  in  Section  1,  the  method  of 
determining  eQ.  In  many  instances  (for  example,  GCG  methods)  this  is  done 

by  an  orthogonal  projection  e0  onto  Sm  with  respect  to  a  suitably  defined 
inner  product.  Thus  if  F  e  Rnxn  is  a  given  SPD  operator,  and  we  define  the 


inner  product 

with  induced  norm  || •  || p .  then 


(v)F  =  C.F-). 


?  =  n  e„. 

o  s„  o 

m 


(6.5) 


where  FT  is  the  orthogonal  projector  of  Rn  onto  S  with  respect  to 
S  m 

m 

(•■•)F. 

To  obtain  a  more  concrete  representation  of  TT  ,  we  assume  that 

m 

Sm  =  Rg(S)  for  some  S  e  Rnxm.  Then  one  easily  has 

n  =  S(StFS)_1  SlF.  (6.6) 

o 

m 

Note  that  by  (6.5)  and  (6,3), 

eQ  =  S(StFS)' 1  sVA'1  r  .  (6.7) 

Hence,  the  prescription  (6.5)  is  practical  only  if  F  is  of  the  form 

F  =  EA  (6.8) 

for  some  E  e  Rnxn.  If  A  is  SPD,  then  an  obvious  and  common  choice  of  E  is 
the  identity  l  (i.e.  F  =  A).  Note  that  in  Section  1  we  also  took  E  =  l. 

In  terms  of  the  element  u  and  ult  equation  (6.5)  is  equivalent  to  the 
condition 


(6.9) 


|u  -  u,|F  <  |w  -  u|F. 

(6.9) 

for  any  element  v  of  the  coset  u0  ♦  Sm.  Again,  if  A  is  SPD  and  F  =  A. 
then  uj  is  also  the  minimizer  over  u0  ♦  Sm  of  the  quadratic  functional 
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If  A  itself  is  SPD,  then  as  previously  observed,  we  can  take  F  =  A.3  In 
this  case  the  contraction  number  of  Theorem  6.1  is 


2TCA-1 )  = 


X2.  (A'1) 
!  _  min 


X^  (A 
max 


-1' 


-1/2 

1  -  — 1 - 1 

= 

k2(A) 

L  2  J 

-|1/2 


where  k2C A)  denotes  the  spectral  number  of  A.  This  should  be  compared  with 
the  contraction  number  (k2(A)  -  1)/(k2(A)  ♦  1)  previously  obtained  by  Chou  and 
Porsching  (1989)  for  this  special  case.  (See  also  Section  1.)  It  is  easy  to  see 
that  2F(A_  1 )  >  (k2(A)  -  l)/(k2(A)  ♦  1).  with  equality  if  and  only  if  A  =  I. 

Next  we  consider  the  case  when  A  is  PR.  Since  A*  is  then  also  PR.  we 
can  take  F  =  A*A.  The  contraction  number  is  2r((A*A)1/2  a-1  (A^A)'1/2).  but 
this  can  be  simplified  by  applying  the  following  lemma  (see  Chou  and  Porsching 
[1990]). 


LEMMA  6.2.  Let  A.P  e  Rnxn  be  respectively  nonsinguiar  and  SPD.  Then 
for  any  real  number  <x, 

jjl  -  oKAtp-lA)1/2  A'1P(Atp-lA)-1/2||2  =  1 1  -  o<p1/2A-1P1/2||2. 

If  A  is  PR.  we  apply  the  lemma  with  P  =  1.  The  result  is  that  the 
contraction  number  ^(A^A)1  /2  A-1  (A*aH /2)  may  be  replaced  by  2T(A_  1 ). 

Finally,  we  turn  to  the  generalized  conjugate  acceleration  procedures  as 
presented  in  Hageman  and  Young  [1981].  Let  Q  e  Rnxn  be  a  nonsingular 
splitting  or  preconditioning  matrix.  Then  (6.1)  can  be  written  as 

(I  -  G)u  =  b.  (6.1 1 ) 

where  G  =  I  -  0_1  A  and  b  =  Q"1  f.  In  the  context  of  the  basic  iterative 
method 

v,  ,  =  Gv^  ♦  b  (v„  =  u  ),  k  =  0.1 .  (6.12) 

k»1  k  00 

the  quantity  r0  =  b  -  (1  -  G)u0  is  known  as  a  "pseudoresiduar.  Note  that 
r0  =  v,  -  v0. 

We  assume  (as  do  Hageman  and  Young  [1981.  p.  341])  that  there  is  a  matrix 
Z  such  that  Z(1  -  G)  is  SPD.  If  Sm  is  the  m  dimensional  (Krylov)  subspace 

spanned  by  the  elements  r0.  (I  -  G)r0 .  and  if  F  =  Z(I  -  G).  then  our 

additive  correction  method  is  equivalent  to  the  generalized  conjugate  gradient 
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method  termed  ORTHOD1R  by  Hageman  and  Young.  If  Z  is  PR.  then  the  nrethod 
becomes  the  so-called  ORTHOMIN  method.  It  follows  immediately  from  the 
definitions  of  Sm  and  F  that  Theorem  6.1  applies  to  ORTHOMIN  with  the 
contraction  number  2f((z(I  -  G))'1/2  z(Z(l  -  G))"1/2). 

The  GCW  generalized  conjugate  gradient  method  (after  Concus  and  Goulb 
[1976]  and  Widlund  [1978])  is  obtained  from  (6.11)  by  choosing  Q  =  1/2(A  ♦  A*) 
under  the  assumption  that  A  is  PR.  In  this  case  we  can  take  F  =  A^O* 1  A  and 
Theorem  6.1  holds  with  respect  to  the  contraction  number 

2f((AtO" 1  A) 1  /2  A' 1  QtA^O"  1  A)'  1  /2).  By  applying  Lemma  6.2  with  P  =  0.  we 
see  that  this  contraction  number  can  be  replaced  by  2f(Q1  /2A-1  Q1  /2). 


7.  Hybrid  Difference  Methods 

In  this  section  we  study  hybrid  difference  methods  for  the  linear  convection 
equation 

uf  ♦  A(x)ux  =  0,  t  >  0.  (7.1) 

subject  to  the  initial  condition, 

u(x,0)  =  ^(x). 

Here  u:  R  -*  Rn.  A:  R  -•  Rn*n,  \p:  R  -»  Rn  are  smooth  functions.  The  matrix 
function  A  is  required  to  be  symmetric,  and  there  exists  a  constant  such 
that  ul  >  A  >  0.  (The  ordering  here  is  that  of  symmetric  matrices.)  By  a 
hybrid  difference  method  we  mean  a  method  that  is  obtained  oy  forming 
weighted  combinations  of  difference  quotients  defining  two  consistent  methods. 

While  there  are  many  ways  to  select  weights  in  the  blending  process,  we  shall 
concentrate  on  a  principle  suggested  by  Porsching  [1989]  for  the  equation  (7.1) 
in  the  case  of  n  =  1.  His  idea  is  to  examine  the  continuous  problem  (7.1)  to  see 
if  a  conservation  of  energy  can  be  found.  If  that  can  be  found,  then  one  tries  to 
create  a  blending  process  so  that  a  similar  energy  conservation  holds  for  the 
corresponding  discrete  case.  Thus  we  start  with  problem  (7.1)  and  derive  a 
conservation  of  weighted  energy  for  it.  Without  loss  of  generality,  we  shall 
assume  the  data  A(x).  \p{x)  given  in  (7.1 )  is  spatially  periodic  with  period  l. 

The  energy  associated  with  (7.1)  is  defined  as 


l2(t)  =  J  (u(x.t),u(x,t))dx. 

0 

where  (•_•)  denotes  the  Euclidean  inner  product.  Note  that  the  energy  so 
defined  is  nothing  but  the  L2-norm  of  the  function  u(-.t)  if  we  define 


1 

<v.w>  s  J  (v.w)dx  for  any  v.w  e  L2(0,4). 


It  is  well  known  that  spatially  periodic  symmetric  hyperbolic  system  has  a 
unique  periodic  solution  (see  Kreiss  (1989,  p.  100]).  Hence  using  the  periodicity 
of  u  and  A.  it  is  easy  to  see  that 

<u.Axu>  =  -2<u,Aux>. 

Thus 

-4-  <u,u>  =  2<u  u> 
dt  t 

=  -2<Au  ,u> 
x 

=  <u.A  u>. 
x 

Hence  I2 ( t)  <  | A  ||  I2(t),  which  implies 

dt  -  11  X 11  oo 


I2(t)  <  expflAxIootJ^CO). 

We  see  that  I2(t)  is  bounded  if  we  are  interested  in  seeking  solution  over  a 
finite  interval  [0 ,T],  In  general,  the  energy  is  not  conserved.  However,  the 
weighted  energy 


M, 

J2(t)  S  J  (u(x.t).A' 1  (x)  u(x.t))dx 


is  conserved.  In  fact 
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J2(t)  =  <u.A_1u>  =  <u.,A_1u>  ♦  <u,A_1u  > 

dt  dt  t  '■ 

=  2<ut,A'1u> 

=  -2<Au  ,A~ 1  U> 
x 

=  0. 

It  is  the  discrete  analog  of  the  weighted  energy  that  we  shall  base  our  weighted 
selection  principle  on.  Let  the  x  -  t  plane  be  covered  by  a  rectangular  mesh 
with  uniform  x  and  t  spacing  h  and  t.  Let  u  be  a  mesh  function  whose 
values  are  in  Rn  and  whose  value  at  the  mesh  (jh.mz)  is  denoted  by  Uj(m). 
When  no  confusion  can  arise  we  shall  suppress  the  dependence  on  j  and/or  m. 
For  any  such  function  we  define  the  common  x-directional  differences: 

S±u  =  u .  . 

x  j*l 

A  v  =  (u  -  u.)/h. 
x  ]*1  J 

V  V  :  (V.  -  V.  ,  )/h. 

X  j  J-1 

A  v  =  (u.  -  u.  ,  )/2h, 

X  j*1  j-1 

S2u  -  ( u .  ,  -  2v.  *  v  ,  )/h2, 

X  ]♦!  j  J*1 

as  well  as  analogous  differences  in  the  t-direction.  The  following  identities 
hold  for  any  mesh  functions  u  and  u: 


V  (u,A  u)  =  (u,s2u)  ♦  (V  u.7  u). 

XX  X  XX 

(7.2) 

V  (u,u) 

X 

=  (u,V  v)  *  (S'u.V  u). 

X  XX 

(7.3) 

2(u.V  u) 

X 

=  V  (u.u)  ♦  h(V  u.V  u). 

X  XX 

(7.4) 

A^u.u) 

=  (u.A^u)  ♦  (A^u.S' u), 

(7.5a) 

1  u. A^u)  = 

A^B^u.u)  -  rCB'1  A(u, Atu). 

(7.5b) 

where  B  =  B(x)  e  Rn*n  is  symmetric. 

These  identities  can  be  proven  first  for  the  case  n  =  1  and  then  for  the 
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general  n  by  applying  the  results  to  each  component.  (7.5b)  follows  from 
(7.5a)  by  setting  u  =  B' 1  v. 

Now  consider  the  following  hybrid  difference  equation  for  (7.1): 

A  v  ♦  AKl  -  9)7  u  *  QA  u]  =  0,  (7.6) 

t  X  X 

where  Aj(m)  =  Aj  =  A(jh)  and  9  is  a  scalar  mesh  function  of  weights.  To  be 
consistent  with  the  continuous  problem,  we  assume  that  v  is  L-periodic  in  j. 
i.e.  u,t|_(m)  =  Uj(m).  Equation  (7.6)  can  be  rewritten  as 

A  v  ♦  A(V  v  ♦  -1-  6h  S2u)  =  0,  (7.7) 

t  x  2  x 

which  clearly  reveals  its  antidiffusive  nature  when  0  >  0.  In  the  next  few 
steps  we  shall  transform  (7.7)  into  a  more  amenable  form  from  which  the 
discrete  energy  can  be  singled  out.  The  idea  is  to  treat  all  new  terms  not 
present  in  the  continuous  case  as  parts  of  the  'source"  term.  By  source  term  we 
mean  the  term  that  tends  to  zero  as  h  and  z  to  zero. 

From  (7.7)  we  have 

(A'V A  u)  *  (A* 1  u,A[V  u  ♦  —  8h  S2u])  =  0. 
t  x  2  x 

Using  (7.2)  to  transform  the  second  term,  we  get 

2(A'Va  u)  ♦  2(u,V  v)  ♦  V  (ehu.A  v)  -  (V  (ehu).v  u)  =  0.  (7.8) 

t  X  X  X  X  X 

Next  we  use  (7.5)  and  (7.4)  to  transform  the  first  two  terms  of  (7.8)  and  the 

identity  V  (0u)  =  (V  e)u  ♦  S*9(V  u)  to  transform  the  fourth  term.  Upon 
x  x  xx 

i  earranging,  we  have 

A  (A*1  u.u)  ♦  V  [(u.u)  ♦  (9hu.A  u)} 
t  x  x 

=  rfA1  A  v,Av)  *  h(V  u,S'e(V  u))  -  h(V  u.V  v)  *  h(V  u.(V  e)v). 

t  t  XXX  XX  XX 

(7.9) 

Adding  h#(Vxu.Vxu)  to  both  sides  of  (7.9).  where  2f  is  a  constant,  we  obtain 
A.(A’1u,v)  ♦  V  [(u.u)  ♦  (Ghu, A  u)]  ♦  h2f(V  u,V  v) 

t  X  XXX 

=  z(L~ 1  A  u,A.  u)  -  h(l  -  Z  -  S'0)(V  u, A  u)  ♦  h(V  G.uV  u).  (7.10) 

t  t  XXX  XX 


Letting 
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b  =  V  0  >  (7.11) 

x 

and 

0  =  z(L']A  u.A  u)  -  h(1  -  V  -  S'u)(V  u.V  u).  (7.12) 

t  t  XXX 

we  see  that 

A  u)  *  V  [(u.u)  ♦  (ehu.A  u)]  ♦  hEf(V  u,V  u)  =  0  +  h(b,uV  u). 

t  x  xxx  x 

(7.13) 


If  we  define  the  weighted  energy  of  the  method  (7.6)  as 

L-1 

P(m)  =  5]  (u.(m).A*1  v  ,(m))h. 

Ti  ]  ]  J 

then  we  can  use  (7.1  3)  to  bound  the  energy 

L-1 

P(m)  =  (u.(m),u.(m))h. 

to  1  1 

provided  that  we  impose  proper  conditions. 

In  the  scalar  case  n  =  1  of  (7.1),  Porsching  [1989]  has  porposed  two 
different  strategies:  global  weight  and  local  weight  selection  procedures. 

However,  we  have  only  been  able  to  generalize  the  global  weight  procedure 
to  higher  dimensional  cases  at  this  moment.  We  describe  how  this  can  be  done. 

Assume  that  the  scalar  weight  function  9  is  independent  of  the  space  idex  j. 

In  this  case  b  =  0.  and  if  we  choose  2f  =  0.  then  (7.13)  reduces  to 

A  u)  ♦  V  [(u.u)  ♦  (0hu,A  u)]  =  0.  (7.14) 

t  x  x 

where 

0  =  z{L']  A  u , A  u)  -  h(1  -  0)(V  u,V  u).  (7.15) 

t  t  xx 

If  we  multiply  (7.1 4)  by  hr,  sum  over  0<j<L-1.  0<k<m,  use  the 
periodicity  of  u.  and  note  that 

y.  A^Z,*1  u.u)hr  =  J2(m  ♦  1)  -  P(0). 
i.k 


we  see  that 
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J2(m  *■  1)  =  3^(0)  *  ^  Q hr.  '  ( 

j.k 

Hence  the  weighted  energy  is  conserved  if 

^  Qht  =  0.  ( 

i.k 

To  enforce  (7.17)  we  observe 

=  (X[V  v  ♦  —  9hS2u],V  u  ♦  3-  9hS2u)  -  ((1  -  0)V  u.V  u ) ,  ( 

h  x  2  x  x  2  x  xx 

where  X  =  t£/h  is  the  mesh  matrix  function  of  "Courant  numbers”.  Thus 


where 


5^  Qhr  =  (D02  +  E0  ♦  F)h2t, 


D  =  Y  (S2u.XS2u). 
4  '  x  x 

i 


E  =  y  (Xh82U.V  u)  ♦  (V  u.V  u] 

™  Y  Y  Y  Y 


F  =  y  (XV  U.V  u)  -  (V  u.V  u). 

XX  XX 


J 

Hence  (7.17)  follows  if  Q  is  a  real  root  of  the  quadratic  function 

q  (0)  =  D02  ♦  EG  ♦  F. 

If  the  "Courant  condition"  p(X)  =  the  spectral  radius  of  xL/h  <  1  holds,  then 

qjO)  =  (XV  u.V  u)  -  (V  u.V  u)  <  0. 

2  *—*  xx  xx  — 

j 

and 


q.(  i )  =  (X1/2  A  u,X1/2A  u)  >  0. 

2  x  x 

j 

Hence  (7.20)  has  a  root  in  the  interval  [0,1]. 

in  the  derivation  so  far,  we  have  assumed  the  scalar  weight  function  0 
depends  on  the  time-spacing,  but  not  on  the  spatial-spacing.  From  now  on,  we 


7.16) 


7.17) 


7.18) 


7.19) 


7.20) 


27 


further  assume  that  0  is  independent  of  the  time  index  k  (i.e.  e  =  constant). 

In  general.  (7.17)  can  no  longer  hold.  However,  if  we  write  (7.18)  in  the  form 

—  =  (XC(1  -  —  e) V  u  ♦  -L  0A  u],[(1  -  -L  0)v  u ♦  —  u])  .  (d  .  q)v  u  y  u) 

h  2  x  2  x  2  x  2  x  x  x 

let  A  =  pr/h,  and  use  2ab  <  E2a2  ♦  e'2b2,  then  we  do  have 

A{|(1  -  e)Vxu)||2  *  2 1 ( 1  -  -L  e)vxu||||l  6Axu||  ♦  [±  eA^upl 

-  (1  -  e)||Vxu|2 

<  A{(  1  -  ^  0)2  (1  ♦  e2) | Vx u || 2  ♦  ©2(i  ♦  e'2)||Vxu||2}  -  (1  -  e)||vxu||2, 

where  |>|  denotes  the  Euclidean  norm.  Thus 

^  Oht  £  q^O)  ^  ||Vxu||2  h2r, 

j  j 

where 

q.(0)  ■  (2  ♦  e  ♦  e~2)02  ♦  [1  -  (1  ♦  e2)A]0  -  (1  ♦  e2)A  -  1.  (7.21) 

Note  that  if  0  <  A  < — - - .  then 

1  ♦  E2 

q,(0)  =  (1  ♦  e2)A  -  1  <  0 
and 

q  (1 )  =  -A-  (2  ♦  e2  ♦  e"2)  >  0. 

1  4 

Thus  there  is  a  root  in  [0,1], 

We  summarize  our  findings  as  follows 

THEOREM  7.1 .  If  A  =  pz/ h  <1/(1  ♦  e2)  and  if  0  is  a  positive  root  of 
the  equation  q]  (0)  =  0  in  (7.21)  for  any  e,  then  0  e  [0.1)  and  J(m)  <  J(0), 

m  >  0. 

THEOREM  7.2.  If  -7-  p(A)  <  1  and  0(m)  is  the  root  of 


lying  in  [0.1],  where 


0(m)32  ♦  E(m)0  ♦  F(m)  =  0. 
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D(m)  =  ^  u  (m).X.(m)S2  u.(m  ). 

4  .  x  j  j  x  j 

J 

E(m)  =  y,  (X.(m)hS2  u.(m),V  u.(m))  +  (V  u.(m),V  u.(m)), 

.  J  x  J  x  J  x  J  *  J 

J 

F(m)  =  2-i  (X.(m)V  u.(m),V  u.Cm))  -  (V  u.(m),V  u.(m)), 

“  J  x  ]  x  J  x  J  x  J 

J 

then  the  weighted  energy  of  (7. 7)  is  conserved,  i.e., 

J(m)  =  J(0)  for  all  m  >  0. 
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